!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!
!
C
*---------------------------------------------------------------------*
c/* Deck CC_CMAT */
*=====================================================================*
      SUBROUTINE CC_CMAT( ICTRAN,  NCTRAN, LISTA, LISTB, LISTC,
     &                    IOPTRES, FILCMA, ICDOTS,CCONS, MXVEC,
     &                    WORK,    LWORK                       )
*---------------------------------------------------------------------*
*
*    Purpose: AO-direct calculation of a linear transformation of three
*             CC amplitude vectors, T^A, T^B and T^C, with the coupled
*             cluster C matrix (second derivative of the CC jacobian 
*             with respect to the amplitudes)
*          
*             The linear transformations are calculated for a list
*             of T^A, T^B and T^C vectors: 
*
*                LISTA       -- type of T^A vectors
*                LISTB       -- type of T^B vectors
*                LISTC       -- type of T^C vectors
*                ICTRAN(1,*) -- indeces of T^A vectors
*                ICTRAN(2,*) -- indeces of T^B vectors
*                ICTRAN(3,*) -- indeces of T^C vectors
*                ICTRAN(4,*) -- indeces or addresses of result vectors
*                NCTRAN      -- number of requested transformations
*                FILCMA      -- file name / list type of result vectors
*                               or list type of vectors to be dotted on
*                ICDOTS      -- indeces of vectors to be dotted on
*                CCONS       -- contains the dot products on return
*
*    return of the result vectors:
*
*           IOPTRES = 0 :  all result vectors are written to a direct
*                          access file, FILCMA is used as file name
*                          the start addresses of the vectors are
*                          returned in ICTRAN(4,*)
*
*           IOPTRES = 1 :  the vectors are kept and returned in WORK
*                          if possible, start addresses returned in
*                          ICTRAN(4,*). N.B.: if WORK is not large
*                          enough iopt is automatically reset to 0!!!
*
*           IOPTRES = 3 :  each result vector is written to its own
*                          file by a call to CC_WRRSP, FILCMA is used
*                          as list type and ICTRAN(4,*) as list index
*                          NOTE that ICTRAN(4,*) is in this case input!
*
*           IOPTRES = 4 :  each result vector is added to a vector on
*                          file by a call to CC_WRRSP, FILCMA is used
*                          as list type and ICTRAN(4,*) as list index
*                          NOTE that ICTRAN(4,*) is in this case input!
*
*           IOPTRES = 5 :  the result vectors are dotted on a array
*                          of vectors, the type of the arrays given
*                          by FILCMA and the indeces from ICDOTS
*                          the result of the dot products is returned
*                          in the CCONS array
*
*
*     Written by Christof Haettig, april-june 1997.
*
*     included CC-R12: Christian Neiss, nov. 2005
*
*=====================================================================*
#if defined (IMPLICIT_NONE)
      IMPLICIT NONE  
#else
#  include "implicit.h"
#endif
#include "priunit.h"
#include "ccsdinp.h"
#include "ccsdsym.h"
#include "maxorb.h"
#include "mxcent.h"
#include "ccorb.h"
#include "cciccset.h"
#include "cbieri.h"
#include "distcl.h"
#include "iratdef.h"
#include "eritap.h"
#include "ccisao.h"
#include "ccfield.h"
#include "aovec.h"
#include "blocks.h"
#include "r12int.h"
#include "ccr12int.h"

* local parameters:
      CHARACTER MSGDBG*(17)
      PARAMETER (MSGDBG='[debug] CC_CMAT> ')

      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.)

      INTEGER KDUM
      PARAMETER( KDUM = +99 999 999 ) ! dummy address for work space
      INTEGER ISYM0
      PARAMETER( ISYM0 = 1 ) ! symmetry of the reference state
      INTEGER ISYOVOV
      PARAMETER( ISYOVOV = 1 ) ! symmetry of (ia|jb) integrals 
      
      INTEGER LUCMAT

      CHARACTER*(*) LISTA, LISTB, LISTC, FILCMA
      INTEGER IOPTRES
      INTEGER NCTRAN, MXVEC, LWORK
      INTEGER ICTRAN(4,NCTRAN)
      INTEGER ICDOTS(MXVEC,NCTRAN)

      DOUBLE PRECISION WORK(LWORK) 
      DOUBLE PRECISION ZERO, ONE, TWO, HALF
      DOUBLE PRECISION CCONS(MXVEC,NCTRAN)
      PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0, HALF = 0.5d0)


      CHARACTER*(10) MODEL, MODELW
      CHARACTER*(3) LISTXA, LISTXB, LISTXC
      CHARACTER*(8) CBAFIL, DBAFIL
      INTEGER INDEXA(MXCORB_CC)
      INTEGER IOPTW, IOPT, ITRAN, IADRTH
      INTEGER LENALL, LEN, IOFFCD, ICYCLE, ILLL, IERROR
      INTEGER KEND0, LEND0, KENDE2, LENDE2, KENDE3, LENDE3, KEND, LEND
      INTEGER KLIAJB, KT2AMPA, KT2AMPC, KT2AMPB, KTHETA1, KTHETA2
      INTEGER KT1AMPA, KT1AMPB, KT1AMPC, KFOCKA, KFOCKB, KFOCKC
      INTEGER KFABVV, KFABOO, KFBCVV, KFBCOO, KFCAVV, KFCAOO, KTHETA0
      INTEGER KX4O, KXIAJB, KKINTC, KCBAR, KDBAR, LUCBAR, LUDBAR
      INTEGER IDLSTA, IDLSTB, IDLSTC, ISYRES, ISYMAB, ISYCBAR, ISYDBAR
      INTEGER ISYMTA, ISYMTB, ISYMTC, ISYMFA, ISYMFB, ISYMFC, ISYX4O
      INTEGER ISYFAB, ISYFBC, ISYFCA, ISYMD1, NTOT, ISYMKC
      INTEGER NTOSYM, KCCFB1, KINDXB, KFREE, LFREE, KENDSV, LWRKSV
      INTEGER KODCL1, KODBC1, KRDBC1, KODPP1, KRDPP1, KRECNR
      INTEGER KODCL2, KODBC2, KRDBC2, KODPP2, KRDPP2, NUMDIS
      INTEGER IDEL2, ISYDEL, IDEL, KXINT
      INTEGER KT1AMP0, KLAMDH0, KLAMDP0, KLAMDHA, KLAMDPA
      INTEGER KLAMDHB, KLAMDPB, KLAMDPC, KLAMDHC
      INTEGER KENDF1, LENDF1, KENDF2, LENDF2, KEND1, LEND1
      INTEGER KEND2, LEND2, KEND3, LEND3, KENDA3, LENDA3
      INTEGER KENDB3, LENDB3, IOPTTCME
      INTEGER IDUM, IAN, KTR12AM, KVINT, LUNIT
      
      DOUBLE PRECISION XNORM 




* external functions:
      INTEGER ILSTSYM

      DOUBLE PRECISION DDOT 
  
      CALL QENTER('CC_CMAT')

*---------------------------------------------------------------------*
* begin:
*---------------------------------------------------------------------*
      IF (LOCDBG) THEN
        Call AROUND('ENTERED CC_CMAT')
        IF (DIRECT) WRITE(LUPRI,'(/1X,A)') 'AO direct transformation'
        WRITE (LUPRI,*) MSGDBG,'LISTA : ',LISTA(1:2)
        WRITE (LUPRI,*) MSGDBG,'LISTB : ',LISTB(1:2)
        WRITE (LUPRI,*) MSGDBG,'LISTC : ',LISTC(1:2)
        WRITE (LUPRI,*) MSGDBG,'FILCMA: ',FILCMA
        WRITE (LUPRI,*) MSGDBG,'NCTRAN: ',NCTRAN
        WRITE (LUPRI,*) MSGDBG,'IOPTRES:',IOPTRES
        WRITE (LUPRI,*) MSGDBG,'MXVEC:',MXVEC
        WRITE (LUPRI,*) MSGDBG,'LWORK:',LWORK
        CALL FLSHFO(LUPRI)
      END IF
      
      IF (CCSDT) THEN
        WRITE(LUPRI,'(/1x,a)') 'C matrix transformations not '
     &          //'implemented for triples yet...'
        CALL QUIT('Triples not implemented for C '//
     &            'matrix transformations')
      END IF

      IF ( .not. (CCS .or. CC2 .or. CCSD) ) THEN
        WRITE(LUPRI,'(/1x,a)') 'CC_CMAT called for a Coupled Cluster '
     &          //'method not implemented in CC_CMAT...'
        CALL QUIT('Unknown CC method in CC_CMAT.')
      END IF

      IF (    LISTA(1:1).NE.'R' 
     &   .OR. LISTB(1:1).NE.'R'
     &   .OR. LISTC(1:1).NE.'R' ) THEN
        WRITE(LUPRI,*) 'LISTA, LISTB and LISTC must refer to',
     &                    ' t-amplituded vectors in CC_CMAT.'
        CALL QUIT('Illegal LISTA or LISTB or LISTC in CC_CMAT.')
      END IF

      IF (ISYMOP .NE. 1) THEN
        WRITE(LUPRI,*) 'ISYMOP = ',ISYMOP
        WRITE(LUPRI,*) 'CC_CMAT is not implemented for ISYMOP.NE.1'
        CALL QUIT('CC_CMAT is not implemented for ISYMOP.NE.1')
      END IF

C     IF (NCTRAN .GT. MAXSIM) THEN
C       WRITE(LUPRI,*) 'NCTRAN = ', NCTRAN
C       WRITE(LUPRI,*) 'MAXSIM = ', MAXSIM
C       WRITE(LUPRI,*) 'number of requested transformation is larger'
C       WRITE(LUPRI,*) 'than the maximum number of allowed ',
C    &                 'simultaneous transformation.'
C       CALL QUIT('Error in CC_CMAT: NCTRAN is larger than MAXSIM.')
C     END IF

* check return option for the result vectors:
      IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN

        LUCMAT = -1
        CALL WOPEN2(LUCMAT,FILCMA,64,0)

      ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN
        CONTINUE
      ELSE IF (IOPTRES.EQ.5) THEN
        IF (MXVEC*NCTRAN.NE.0) CALL DZERO(CCONS,MXVEC*NCTRAN)  
      ELSE
        CALL QUIT('Illegal value of IOPTRES in CC_CMAT.')
      END IF

* construct 'MODEL' string for CC_WRRSP routine and set write option:
      IF (CCS) THEN
         MODELW = 'CCS       '
         IOPTW  = 1
      ELSE IF (CC2) THEN
         MODELW = 'CC2       '
         IOPTW  = 3
      ELSE IF (CCSD) THEN
         MODELW = 'CCSD      '
         IOPTW  = 3
      ELSE
         CALL QUIT('Unknown coupled cluster model in CC_CMAT.')
      END IF



* start of scratch space for the following:
      KEND0 = 1
      LEND0 = LWORK

*=====================================================================*
* Calculate J term and N^5 terms E1 and E2 in a loop over all required
* C matrix transformations
*
* the N^5 terms H, G, I vanish for the C matrix 
*
* All N^5 terms drop out for CCS, the E1, E2 vanish also for CC2.
*=====================================================================*

* memory allocation:
      KLIAJB = KEND0
      KENDE2 = KLIAJB + NT2AM(ISYOVOV)
      LENDE2 = LWORK - KENDE2

      IF (LENDE2.LT.0) THEN
        CALL QUIT('Insufficient memory in CC_CMAT.(1)')
      END IF

* read packed (ia|jb) integrals and calculate L(ia|jb) in place: 
      Call CCG_RDIAJB(WORK(KLIAJB),NT2AM(ISYOVOV))
      
      IOPTTCME = 1
      Call CCSD_TCMEPK(WORK(KLIAJB),ONE,ISYOVOV,IOPTTCME)

*---------------------------------------------------------------------*
* start loop over all C matrix transformations
*---------------------------------------------------------------------*
      IADRTH = 1

      DO ITRAN = 1, NCTRAN
        IDLSTA = ICTRAN(1,ITRAN)
        IDLSTB = ICTRAN(2,ITRAN)
        IDLSTC = ICTRAN(3,ITRAN)

        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG, 'ITRAN :',ITRAN
          WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),IDLSTA
          WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),IDLSTB
          WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),IDLSTC
          CALL FLSHFO(LUPRI)
        END IF
*---------------------------------------------------------------------*
* set symmetries and allocate memory:
*---------------------------------------------------------------------*
        ISYMTA = ILSTSYM(LISTA,IDLSTA)
        ISYMTB = ILSTSYM(LISTB,IDLSTB)
        ISYMTC = ILSTSYM(LISTC,IDLSTC)

        ISYMFA = MULD2H(ISYOVOV,ISYMTA)
        ISYMFB = MULD2H(ISYOVOV,ISYMTB)
        ISYMFC = MULD2H(ISYOVOV,ISYMTC)

        ISYFAB = MULD2H(ISYMFA,ISYMTB)
        ISYFBC = MULD2H(ISYMFB,ISYMTC)
        ISYFCA = MULD2H(ISYMFC,ISYMTA)

        ISYRES = MULD2H(ISYFAB,ISYMTC)

* allocate memory for double excitation response vectors:
* (overlaps with memory for two-electron integrals)
        KENDE3  = KLIAJB + NT2AM(ISYOVOV)
        IF (.NOT.CCS) THEN
          KT2AMPA = KLIAJB
          KT2AMPB = KLIAJB
          KT2AMPC = KLIAJB

          KENDE3  = MAX(KENDE3,KLIAJB + NT2SQ(ISYMTA))
          KENDE3  = MAX(KENDE3,KLIAJB + NT2SQ(ISYMTB))
          KENDE3  = MAX(KENDE3,KLIAJB + NT2SQ(ISYMTC))
        END IF

* allocate memory for the result vector:
        KTHETA1 = KENDE3
        KENDE3  = KTHETA1 + NT1AM(ISYRES)

        IF (.NOT.CCS) THEN
          KTHETA2 = KENDE3
          KENDE3  = KTHETA2 + NT2AM(ISYRES)
        END IF

        KT1AMPA = KENDE3
        KT1AMPB = KT1AMPA + NT1AM(ISYMTA)
        KT1AMPC = KT1AMPB + NT1AM(ISYMTB)
        KFOCKA  = KT1AMPC + NT1AM(ISYMTC)
        KFOCKB  = KFOCKA  + NT1AM(ISYMFA)
        KFOCKC  = KFOCKB  + NT1AM(ISYMFB)
        KFABVV  = KFOCKC  + NT1AM(ISYMFC)
        KFABOO  = KFABVV  + NMATAB(ISYFAB)
        KFBCVV  = KFABOO  + NMATIJ(ISYFAB)
        KFBCOO  = KFBCVV  + NMATAB(ISYFBC)
        KFCAVV  = KFBCOO  + NMATIJ(ISYFBC)
        KFCAOO  = KFCAVV  + NMATAB(ISYFCA)
        KENDE3  = KFCAOO  + NMATIJ(ISYFCA)
        LENDE3  = LWORK - KENDE3

        IF (LENDE3 .LT. 0) THEN
          CALL QUIT('Insufficient work space in CC_CMAT.(2)')
        END IF

*---------------------------------------------------------------------*
* read single excitation part of the response vectors:
*---------------------------------------------------------------------*
* read A response amplitudes:
        IOPT = 1
        Call CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL,
     &                WORK(KT1AMPA),WORK(KDUM))

* read B response amplitudes:
        IOPT = 1
        Call CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL,
     &                WORK(KT1AMPB),WORK(KDUM))

* read C response amplitudes:
        IOPT = 1
        Call CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL,
     &                WORK(KT1AMPC),WORK(KDUM))

*---------------------------------------------------------------------*
* calculate occupied/virtual part of first order Fock matrices:
*---------------------------------------------------------------------*
* calculate tF^A_ia:
        Call CCG_LXD(WORK(KFOCKA), ISYMFA, WORK(KT1AMPA), ISYMTA, 
     &               WORK(KLIAJB), ISYOVOV, 0)
      
* calculate tF^B_ia:
        Call CCG_LXD(WORK(KFOCKB), ISYMFB, WORK(KT1AMPB), ISYMTB,
     &               WORK(KLIAJB), ISYOVOV, 0)

* calculate tF^C_ia:
        Call CCG_LXD(WORK(KFOCKC), ISYMFC, WORK(KT1AMPC), ISYMTC, 
     &               WORK(KLIAJB), ISYOVOV, 0)
      
*---------------------------------------------------------------------*
* calculate double one-index transformed Fock matrices:
*---------------------------------------------------------------------*
* calculate occ/occ block for FAB:
        Call CCG_1ITROO(WORK(KFABOO),ISYFAB,
     &                  WORK(KFOCKB),ISYMFB,
     &                  WORK(KT1AMPA),ISYMTA    )

        IF (LENDE3.LT.NMATIJ(ISYFAB)) 
     &       CALL QUIT('Out of memory in CC_CMAT.')

        Call CCG_1ITROO(WORK(KENDE3),ISYFAB,
     &                  WORK(KFOCKA),ISYMFA,
     &                  WORK(KT1AMPB),ISYMTB    )

        Call DAXPY(NMATIJ(ISYFAB),ONE,WORK(KENDE3),1,WORK(KFABOO),1)

* calculate vir/vir block for FAB:
        Call CCG_1ITRVV(WORK(KFABVV),ISYFAB,
     &                  WORK(KFOCKB),ISYMFB,
     &                  WORK(KT1AMPA),ISYMTA  )

        IF (LENDE3.LT.NMATAB(ISYFAB)) 
     &       CALL QUIT('Out of memory in CC_CMAT.')

        Call CCG_1ITRVV(WORK(KENDE3),ISYFAB,
     &                  WORK(KFOCKA),ISYMFA,
     &                  WORK(KT1AMPB),ISYMTB  )

        Call DAXPY(NMATAB(ISYFAB),ONE,WORK(KENDE3),1,WORK(KFABVV),1)

* calculate occ/occ block for FBC:
        Call CCG_1ITROO(WORK(KFBCOO),ISYFBC,
     &                  WORK(KFOCKB),ISYMFB,
     &                  WORK(KT1AMPC),ISYMTC    )

        IF (LENDE3.LT.NMATIJ(ISYFBC)) 
     &       CALL QUIT('Out of memory in CC_CMAT.')

        Call CCG_1ITROO(WORK(KENDE3),ISYFBC,
     &                  WORK(KFOCKC),ISYMFC,
     &                  WORK(KT1AMPB),ISYMTB    )

        Call DAXPY(NMATIJ(ISYFBC),ONE,WORK(KENDE3),1,WORK(KFBCOO),1)

* calculate vir/vir block for FBC:
        Call CCG_1ITRVV(WORK(KFBCVV),ISYFBC,
     &                  WORK(KFOCKB),ISYMFB,
     &                  WORK(KT1AMPC),ISYMTC  )

        IF (LENDE3.LT.NMATAB(ISYFBC)) 
     &       CALL QUIT('Out of memory in CC_CMAT.')

        Call CCG_1ITRVV(WORK(KENDE3),ISYFBC,
     &                  WORK(KFOCKC),ISYMFC,
     &                  WORK(KT1AMPB),ISYMTB  )

        Call DAXPY(NMATAB(ISYFBC),ONE,WORK(KENDE3),1,WORK(KFBCVV),1)

* calculate occ/occ block for FCA:
        Call CCG_1ITROO(WORK(KFCAOO),ISYFCA,
     &                  WORK(KFOCKC),ISYMFC,
     &                  WORK(KT1AMPA),ISYMTA    )

        IF (LENDE3.LT.NMATIJ(ISYFCA))
     &       CALL QUIT('Out of memory in CC_CMAT.')

        Call CCG_1ITROO(WORK(KENDE3),ISYFCA,
     &                  WORK(KFOCKA),ISYMFA,
     &                  WORK(KT1AMPC),ISYMTC    )

        Call DAXPY(NMATIJ(ISYFCA),ONE,WORK(KENDE3),1,WORK(KFCAOO),1)

* calculate vir/vir block for FCA:
        Call CCG_1ITRVV(WORK(KFCAVV),ISYFCA,
     &                  WORK(KFOCKC),ISYMFC,
     &                  WORK(KT1AMPA),ISYMTA  )

        IF (LENDE3.LT.NMATAB(ISYFCA))
     &       CALL QUIT('Out of memory in CC_CMAT.')

        Call CCG_1ITRVV(WORK(KENDE3),ISYFCA,
     &                  WORK(KFOCKA),ISYMFA,
     &                  WORK(KT1AMPC),ISYMTC  )

        Call DAXPY(NMATAB(ISYFCA),ONE,WORK(KENDE3),1,WORK(KFCAVV),1)

*---------------------------------------------------------------------*
* J term
*---------------------------------------------------------------------*
        IF (LENDE3.LT.NT1AM(ISYRES))
     &       CALL QUIT('Out of memory in CC_CMAT.')

* initialize single-excitation part of the result vector:
        CALL DZERO(WORK(KTHETA1),NT1AM(ISYRES))

* calculate triple one-index transformed contribution T^C x F^AB:
        IOPT  = 1
        Call CCG_1ITRVO( WORK(KENDE3), ISYRES,
     &                   WORK(KFABOO), WORK(KFABVV), ISYFAB,
     &                   WORK(KT1AMPC), ISYMTC, IOPT          )

        CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KENDE3),1,WORK(KTHETA1),1)


* calculate triple one-index transformed contribution T^A x F^BC:
        IOPT  = 1
        Call CCG_1ITRVO( WORK(KENDE3), ISYRES,
     &                   WORK(KFBCOO), WORK(KFBCVV), ISYFBC,
     &                   WORK(KT1AMPA), ISYMTA, IOPT          )

        CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KENDE3),1,WORK(KTHETA1),1)


* calculate triple one-index transformed contribution T^B x F^CA:
        IOPT  = 1
        Call CCG_1ITRVO( WORK(KENDE3), ISYRES,
     &                   WORK(KFCAOO), WORK(KFCAVV), ISYFCA,
     &                   WORK(KT1AMPB), ISYMTB, IOPT          )

        CALL DAXPY(NT1AM(ISYRES),HALF,WORK(KENDE3),1,WORK(KTHETA1),1)

        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG,'THETA1 after J term:'
          CALL CC_PRP(WORK(KTHETA1),WORK(KDUM),ISYRES,1,0)
          CALL FLSHFO(LUPRI)
        END IF

*---------------------------------------------------------------------*
* initialize double-excitation part of the result vector:
*---------------------------------------------------------------------*
        IF (.NOT.CCS) THEN
          CALL DZERO(WORK(KTHETA2),NT2AM(ISYRES))
        END IF
 
*---------------------------------------------------------------------*
* E1 and E2 terms 
*---------------------------------------------------------------------*
        IF (.NOT. (CCS.OR.CC2.OR.CCSTST) ) THEN

* check available scratch space:
          IF (     LENDE3 .LT. NT2AM(ISYMTA)
     &        .OR. LENDE3 .LT. NT2AM(ISYMTB) 
     &        .OR. LENDE3 .LT. NT2AM(ISYMTC) ) THEN
            CALL QUIT('Insufficient work space in CC_BMAT.(3)')
          END IF

* read double excitation part of T^C vector and unpack to full matrix:
          IOPT = 2
          CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL,
     &                  WORK(KDUM),WORK(KENDE3))
          CALL CCLR_DIASCL(WORK(KENDE3),TWO,ISYMTC)
          CALL CC_T2SQ(WORK(KENDE3),WORK(KT2AMPC),ISYMTC)

* calculate the contribution to THETA2:
          CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPC),WORK(KFABVV),
     &                 WORK(KFABOO),WORK(KENDE3),LENDE3,ISYMTC,ISYFAB)


* read double excitation part of T^A vector and unpack to full matrix:
          IOPT = 2
          CALL CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL,
     &                  WORK(KDUM),WORK(KENDE3))
          CALL CCLR_DIASCL(WORK(KENDE3),TWO,ISYMTA)
          CALL CC_T2SQ(WORK(KENDE3),WORK(KT2AMPA),ISYMTA)

* calculate the contribution to THETA2:
          CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPA),WORK(KFBCVV),
     &                 WORK(KFBCOO),WORK(KENDE3),LENDE3,ISYMTA,ISYFBC)


* read double excitation part of T^B vector and unpack to full matrix:
          IOPT = 2
          CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL,
     &                  WORK(KDUM),WORK(KENDE3))
          CALL CCLR_DIASCL(WORK(KENDE3),TWO,ISYMTB)
          CALL CC_T2SQ(WORK(KENDE3),WORK(KT2AMPB),ISYMTB)

* calculate the contribution to THETA2:
          CALL CCRHS_E(WORK(KTHETA2),WORK(KT2AMPB),WORK(KFCAVV),
     &                 WORK(KFCAOO),WORK(KENDE3),LENDE3,ISYMTB,ISYFCA)


* restore L(ia|jb) integrals: 
          IF (ITRAN.LT.NCTRAN) THEN
            Call CCG_RDIAJB(WORK(KLIAJB),NT2AM(ISYOVOV))
            IOPTTCME = 1
            Call CCSD_TCMEPK(WORK(KLIAJB),ONE,ISYOVOV,IOPTTCME)
          END IF

        END IF

*---------------------------------------------------------------------*
* save result vector:
*---------------------------------------------------------------------*
      IF (LOCDBG) THEN
        WRITE (LUPRI,*) MSGDBG,'THETA after E term:'
        CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1)
        WRITE (LUPRI,*) MSGDBG,'IOPTRES:',IOPTRES
        WRITE (LUPRI,*) MSGDBG,'FILCMA:',FILCMA
        WRITE (LUPRI,*) MSGDBG,'LUCMAT:',LUCMAT
        WRITE (LUPRI,*) MSGDBG,'ICTRAN(4,ITRAN):',ICTRAN(4,ITRAN)
        CALL FLSHFO(LUPRI)
      END IF

      KTHETA0 = -999999


      IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN

*       write to a common direct acces file,
*       store start address in ICTRAN(4,ITRAN)

        ICTRAN(4,ITRAN) = IADRTH

        CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA1),IADRTH,NT1AM(ISYRES))
        IADRTH = IADRTH + NT1AM(ISYRES)

        IF (.NOT.CCS) THEN
          CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA2),IADRTH,NT2AM(ISYRES))
          IADRTH = IADRTH + NT2AM(ISYRES)
        END IF


      ELSE IF (IOPTRES .EQ. 3) THEN

*       write to a sequential file by call to CC_WRRSP,
*       use FILCMA as LIST type and ICTRAN(4,ITRAN) as index

        CALL CC_WRRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW,
     &                WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
     &                WORK(KENDE3),LENDE3)

      ELSE IF (IOPTRES .EQ. 4) THEN

*       add to a vector on a sequential file by call to CC_WARSP,
*       use FILCMA as LIST type and ICTRAN(4,ITRAN) as index

        CALL CC_WARSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW,
     &                WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
     &                WORK(KENDE3),LENDE3)

      ELSE IF (IOPTRES .EQ. 5) THEN
       
*       calculate required dot products

        CALL CCDOTRSP(ICDOTS,CCONS,IOPTW,FILCMA,ITRAN,NCTRAN,MXVEC,
     &                WORK(KTHETA1),WORK(KTHETA2),ISYRES,
     &                WORK(KENDE3),LENDE3)

      ELSE
        CALL QUIT('illegal value for IOPT in CC_CMAT.')
      END IF

* end of loop over C matrix transformations:
      END DO

*=====================================================================*
* end of J, E1 and E2 terms. 
*=====================================================================*

      IF (LOCDBG) THEN
        WRITE (LUPRI,*) MSGDBG,'J and E (CC2/CCSD) term finished...'
        CALL FLSHFO(LUPRI)
      END IF

* that's it for CCS:
      IF (CCS.OR.CCSTST) GOTO 8888

*=====================================================================*
* F term: requires AO integrals...
*
*         Loop over AO integral shells 
*            Loop over C matrix transformations 
*               Loop over AO integral distributions
*  
*                  Calculation of F term contributions
*
*               End loop
*            End loop
*         End loop
*
*  F term drops out for CCS.
*
*=====================================================================*
      IF (.NOT. (CCS.OR.CCSTST) ) THEN


*---------------------------------------------------------------------*
* initialize integral calculation
*---------------------------------------------------------------------*

        KEND = KEND0
        LEND = LEND0

        IF (DIRECT) THEN
           NTOSYM = 1

           IF (HERDIR) THEN
             CALL HERDI1(WORK(KEND),LEND,IPRERI)
           ELSE
             KCCFB1 = KEND
             KINDXB = KCCFB1 + MXPRIM*MXCONT
             KEND   = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT
             LEND   = LWORK  - KEND
             CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2,
     *                   KODPP1,KODPP2,KRDPP1,KRDPP2,
     *                   KFREE,LFREE,KEND,WORK(KCCFB1),WORK(KINDXB),
     *                   WORK(KEND),LEND,IPRERI)
 
             KEND = KFREE
             LEND = LFREE
           END IF

           KENDSV = KEND
           LWRKSV = LEND
        ELSE
           NTOSYM = NSYM
        END IF

*---------------------------------------------------------------------*
* start loop over AO integrals shells:
*---------------------------------------------------------------------*
      DO ISYMD1 = 1, NTOSYM

        IF (DIRECT) THEN
          IF (HERDIR) THEN
             NTOT = MAXSHL
          ELSE
             NTOT = MXCALL
          ENDIF                   
        ELSE
          NTOT = NBAS(ISYMD1)
        END IF

        DO ILLL = 1, NTOT

          IF (DIRECT) THEN
            KEND = KENDSV
            LEND = LWRKSV

            IF (HERDIR) THEN
               CALL HERDI2(WORK(KEND),LEND,INDEXA,ILLL,NUMDIS, 
     &                     IPRINT) 
            ELSE
              CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0,
     *                    WORK(KODCL1),WORK(KODCL2),
     *                    WORK(KODBC1),WORK(KODBC2),
     *                    WORK(KRDBC1),WORK(KRDBC2),
     *                    WORK(KODPP1),WORK(KODPP2),
     *                    WORK(KRDPP1),WORK(KRDPP2),
     *                    WORK(KCCFB1),WORK(KINDXB),
     *                    WORK(KEND), LEND,IPRERI)
            END IF

            KRECNR = KEND
            KEND   = KRECNR + (NBUFX(0) - 1)/IRAT + 1
            LEND   = LWORK - KEND
 
            IF (LEND .LT. 0) THEN
              CALL QUIT('Insufficient work space in CC_CMAT. (1a)')
            END IF

          ELSE
            NUMDIS = 1
          END IF


*---------------------------------------------------------------------*
*         start loop over C matrix transformations:
*---------------------------------------------------------------------*
          DO ITRAN = 1, NCTRAN

            IDLSTA = ICTRAN(1,ITRAN)
            IDLSTB = ICTRAN(2,ITRAN)
            IDLSTC = ICTRAN(3,ITRAN)

            ISYMTA = ILSTSYM(LISTA,IDLSTA)
            ISYMTB = ILSTSYM(LISTB,IDLSTB)
            ISYMTC = ILSTSYM(LISTC,IDLSTC)

            ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),ISYMTC)

* single excitation parts of coupled cluster vectors:
            KT1AMP0 = KEND
            KT1AMPA = KT1AMP0 + NT1AM(ISYM0)
            KT1AMPB = KT1AMPA + NT1AM(ISYMTA)
            KT1AMPC = KT1AMPB + NT1AM(ISYMTB)
            KENDF1  = KT1AMPC + NT1AM(ISYMTC)

* Lambda-hole matrices:
            KLAMDH0 = KENDF1
            KLAMDHA = KLAMDH0 + NGLMDT(ISYM0)
            KLAMDHB = KLAMDHA + NGLMDT(ISYMTA)
            KLAMDHC = KLAMDHB + NGLMDT(ISYMTB)
            KENDF1  = KLAMDHC + NGLMDT(ISYMTC)

* Lambda-particle matrices:
            KLAMDP0 = KENDF1
            KLAMDPA = KLAMDP0 + NGLMDT(ISYM0)
            KLAMDPB = KLAMDPA + NGLMDT(ISYMTA)
            KLAMDPC = KLAMDPB + NGLMDT(ISYMTB)
            KENDF1  = KLAMDPC + NGLMDT(ISYMTC)

* the result vector:
            KTHETA1 = KENDF1
            KTHETA2 = KTHETA1 + NT1AM(ISYRES)
            KENDF1  = KTHETA2 + NT2AM(ISYRES)
            LENDF1  = LWORK - KENDF1

            IF (LENDF1 .LT. 0) THEN
              CALL QUIT('Insufficient memory in CC_CMAT.(1F)')
            END IF

  
* read coupled cluster vectors:
            IOPT = 1
            CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,
     &                    WORK(KT1AMP0),WORK(KDUM))

            IOPT = 1
            CALL CC_RDRSP(LISTA,IDLSTA,ISYMTA,IOPT,MODEL,
     &                    WORK(KT1AMPA),WORK(KDUM))

            IOPT = 1
            CALL CC_RDRSP(LISTB,IDLSTB,ISYMTB,IOPT,MODEL,
     &                    WORK(KT1AMPB),WORK(KDUM))

            IOPT = 1
            CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL,
     &                    WORK(KT1AMPC),WORK(KDUM))

* calculate unperturbed Lambda matrices:
            CALL LAMMAT(WORK(KLAMDP0),WORK(KLAMDH0),WORK(KT1AMP0),
     &                  WORK(KENDF1),LENDF1)

* calculate response Lambda matrices:
            CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPA),WORK(KLAMDH0),
     &                       WORK(KLAMDHA),WORK(KT1AMPA),ISYMTA)

            CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPB),WORK(KLAMDH0),
     &                       WORK(KLAMDHB),WORK(KT1AMPB),ISYMTB)

            CALL CCLR_LAMTRA(WORK(KLAMDP0),WORK(KLAMDPC),WORK(KLAMDH0),
     &                       WORK(KLAMDHC),WORK(KT1AMPC),ISYMTC)


* read result vector:
            IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN
              CALL GETWA2(LUCMAT,FILCMA,WORK(KTHETA1),
     &                    ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES))
            ELSE IF (IOPTRES.EQ.3) THEN
              IOPT = 3
              CALL CC_RDRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPT,MODEL,
     &                      WORK(KTHETA1),WORK(KTHETA2))
            ELSE IF (IOPTRES.EQ.4) THEN
              CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) )
              CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) )
            ELSE IF (IOPTRES.EQ.5) THEN
              CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) )
              CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) )
            ELSE
              CALL QUIT('Error in CC_CMAT.')
            END IF

*---------------------------------------------------------------------*
*        loop over number of distributions on the disk:
*---------------------------------------------------------------------*
          DO IDEL2  = 1, NUMDIS

            IF (DIRECT) THEN
              IDEL   = INDEXA(IDEL2)
              IF (NOAUXB) THEN
                IDUM = 1
                CALL IJKAUX(IDEL,IDUM,IDUM,IDUM)
              END IF
              ISYDEL = ISAO(IDEL)
            ELSE
              IDEL   = IBAS(ISYMD1) + ILLL
              ISYDEL = ISYMD1
            END IF

*           read AO integral distribution and calculate integrals with
*           one index transformed to occupied MO (particle):

            KXINT  = KENDF1
            KENDF2 = KXINT  + NDISAO(ISYDEL)
            LENDF2 = LWORK - KENDF2

            IF (LENDF2 .LT. 0) THEN
              CALL QUIT('Insufficient work space in CC_CMAT. (3)')
            END IF

            CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KENDF2),LENDF2,
     &                  WORK(KRECNR),DIRECT)
*.....................................................................*

* set option for CC_MOFCON routine:
            IOPT = 3

*.....................................................................*
            CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2),
     &                      WORK(KLAMDPA),WORK(KLAMDHA),
     &                      WORK(KLAMDPB),WORK(KLAMDHB),
     &                      WORK(KLAMDPC),WORK(KLAMDH0),
     &                      ISYMTA,ISYMTB,ISYMTC,ISYM0,
     &                      WORK(KENDF2),LENDF2,IDEL,
     &                      ISYDEL,ISYRES,ISYM0,IOPT)

            CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2),
     &                      WORK(KLAMDPA),WORK(KLAMDHA),
     &                      WORK(KLAMDPB),WORK(KLAMDHB),
     &                      WORK(KLAMDP0),WORK(KLAMDHC),
     &                      ISYMTA,ISYMTB,ISYM0,ISYMTC,
     &                      WORK(KENDF2),LENDF2,IDEL,
     &                      ISYDEL,ISYRES,ISYM0,IOPT)

*.....................................................................*
            CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2),
     &                      WORK(KLAMDPB),WORK(KLAMDHB),
     &                      WORK(KLAMDPC),WORK(KLAMDHC),
     &                      WORK(KLAMDP0),WORK(KLAMDHA),
     &                      ISYMTB,ISYMTC,ISYM0,ISYMTA,
     &                      WORK(KENDF2),LENDF2,IDEL,
     &                      ISYDEL,ISYRES,ISYM0,IOPT)

            CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2),
     &                      WORK(KLAMDPB),WORK(KLAMDHB),
     &                      WORK(KLAMDP0),WORK(KLAMDH0),
     &                      WORK(KLAMDPC),WORK(KLAMDHA),
     &                      ISYMTB,ISYM0,ISYMTC,ISYMTA,
     &                      WORK(KENDF2),LENDF2,IDEL,
     &                      ISYDEL,ISYRES,ISYM0,IOPT)

*.....................................................................*
            CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2),
     &                      WORK(KLAMDPB),WORK(KLAMDHB),
     &                      WORK(KLAMDPC),WORK(KLAMDHC),
     &                      WORK(KLAMDPA),WORK(KLAMDH0),
     &                      ISYMTB,ISYMTC,ISYMTA,ISYM0,
     &                      WORK(KENDF2),LENDF2,IDEL,
     &                      ISYDEL,ISYRES,ISYM0,IOPT)
 
            CALL CC_MOFCON2(WORK(KXINT),WORK(KTHETA2),
     &                      WORK(KLAMDPB),WORK(KLAMDHB),
     &                      WORK(KLAMDP0),WORK(KLAMDH0),
     &                      WORK(KLAMDPA),WORK(KLAMDHC),
     &                      ISYMTB,ISYM0,ISYMTA,ISYMTC,
     &                      WORK(KENDF2),LENDF2,IDEL,
     &                      ISYDEL,ISYRES,ISYM0,IOPT)
 
*.....................................................................*


          END DO ! IDEL2
*---------------------------------------------------------------------*
*         end of the loop over integral distributions:
*---------------------------------------------------------------------*
          IF (LOCDBG) THEN
            WRITE (LUPRI,*) MSGDBG,'THETA after F term:'
            WRITE (LUPRI,*) MSGDBG, 'ITRAN :',ITRAN
            WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),
     &           ICTRAN(1,ITRAN)
            WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),
     &           ICTRAN(2,ITRAN)
            WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),
     &           ICTRAN(3,ITRAN)
            CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1)
            WRITE (LUPRI,*) MSGDBG,'IOPTRES:',IOPTRES
            WRITE (LUPRI,*) MSGDBG,'FILCMA:',FILCMA
            WRITE (LUPRI,*) MSGDBG,'LUCMAT:',LUCMAT
            WRITE (LUPRI,*) MSGDBG,'ICTRAN(4,ITRAN):',ICTRAN(4,ITRAN)
            CALL FLSHFO(LUPRI)
          END IF

          KTHETA0 = -999999

          IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN
            CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA1),
     &                  ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES))
c         ELSE IF (IOPTRES.EQ.3) THEN
c           CALL CC_WRRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW,
c    &                    WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
c    &                    WORK(KENDF1),LENDF1)
          ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN
            CALL CC_WARSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW,
     &                    WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
     &                    WORK(KENDF1),LENDF1)
          ELSE IF (IOPTRES.EQ.5) THEN
            CALL CCDOTRSP(ICDOTS,CCONS,IOPTW,FILCMA,ITRAN,NCTRAN,MXVEC,
     &                    WORK(KTHETA1),WORK(KTHETA2),ISYRES,
     &                    WORK(KENDF1),LENDF1)
          ELSE
            CALL QUIT('Error in CC_CMAT.')
          END IF

*---------------------------------------------------------------------*
*         end of the loop over C matrix transformations:
*---------------------------------------------------------------------*
          END DO ! ITRAN
       END DO ! ILLL
      END DO ! ISYMD1
*=====================================================================*
* End of Loop over AO-integrals
*=====================================================================*

        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG,'F term section finished...'
          CALL FLSHFO(LUPRI)
        END IF

      END IF ! (.NOT.CCS)
*=====================================================================*
* end of F term section
*=====================================================================*

* that's it for CC2 
      IF (CC2) GOTO 8888


*=====================================================================*
* Calculate the N^6 terms A, B, C and D in a loop over all required
* C matrix transformations.
*
* All N^6 terms drop out for CCS and CC2.
*=====================================================================*
      IF (.NOT. (CCS .OR. CC2 .OR. CCSTST) ) THEN


*---------------------------------------------------------------------*
* loop over all C matrix transformations:
*---------------------------------------------------------------------*
      DO ITRAN = 1, NCTRAN
        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG, 'N^6 term section:'
          WRITE (LUPRI,*) MSGDBG, 'ITRAN :',ITRAN
          WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),
     &         ICTRAN(1,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),
     &         ICTRAN(2,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),
     &         ICTRAN(3,ITRAN)
          CALL FLSHFO(LUPRI)
        END IF

        ISYMTA = ILSTSYM(LISTA,ICTRAN(1,ITRAN))
        ISYMTB = ILSTSYM(LISTC,ICTRAN(2,ITRAN))
        ISYMTC = ILSTSYM(LISTB,ICTRAN(3,ITRAN))
        ISYRES = MULD2H(MULD2H(ISYMTA,ISYMTB),MULD2H(ISYMTC,ISYOVOV))

* read result vector:
        KTHETA1 = KEND0
        KTHETA2 = KTHETA1 + NT1AM(ISYRES)
        KEND1   = KTHETA2 + NT2AM(ISYRES)
        LEND1   = LWORK - KEND1

        IF (LEND1.LT.0) THEN
          CALL QUIT('Insufficient memory in CC_CMAT. (1b)')
        END IF

        IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN
          CALL GETWA2(LUCMAT,FILCMA,WORK(KTHETA1),
     &                ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES))
        ELSE IF (IOPTRES.EQ.3) THEN
          IOPT = 3
          CALL CC_RDRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPT,MODEL,
     &                  WORK(KTHETA1),WORK(KTHETA2))
        ELSE IF (IOPTRES.EQ.4) THEN
          CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) )
          CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) )
        ELSE IF (IOPTRES.EQ.5) THEN
          CALL DZERO( WORK(KTHETA1), NT1AM(ISYRES) )
          CALL DZERO( WORK(KTHETA2), NT2AM(ISYRES) )
        ELSE
          CALL QUIT('Error in CC_CMAT.')
        END IF

* loop over cyclic permutations of A, B and C:
        DO ICYCLE = 1, 3
          IF (ICYCLE.EQ.1) THEN
             IDLSTA = ICTRAN(1,ITRAN)
             IDLSTB = ICTRAN(2,ITRAN)
             IDLSTC = ICTRAN(3,ITRAN)
           
             LISTXA = LISTA
             LISTXB = LISTB
             LISTXC = LISTC
          ELSE IF (ICYCLE.EQ.2) THEN
             IDLSTB = ICTRAN(1,ITRAN)
             IDLSTC = ICTRAN(2,ITRAN)
             IDLSTA = ICTRAN(3,ITRAN)
             
             LISTXB = LISTA
             LISTXC = LISTB
             LISTXA = LISTC
          ELSE IF (ICYCLE.EQ.3) THEN
             IDLSTC = ICTRAN(1,ITRAN)
             IDLSTA = ICTRAN(2,ITRAN)
             IDLSTB = ICTRAN(3,ITRAN)
             
             LISTXC = LISTA
             LISTXA = LISTB
             LISTXB = LISTC
          ELSE
             CALL QUIT('Error in CC_CMAT.')
          END IF

          ISYMTA = ILSTSYM(LISTXA,IDLSTA)
          ISYMTB = ILSTSYM(LISTXB,IDLSTB)
          ISYMTC = ILSTSYM(LISTXC,IDLSTC)
          ISYMAB = MULD2H(ISYMTA,ISYMTB)


*---------------------------------------------------------------------*
* read single excitation parts of T^A and T^B and keep them in 
* memory during the loop:
*---------------------------------------------------------------------*
        KT1AMPA  = KEND1
        KT1AMPB  = KT1AMPA + NT1AM(ISYMTA)
        KEND2    = KT1AMPB + NT1AM(ISYMTB)
        LEND2    = LWORK - KEND2

        IF (LEND2 .LE. 0) THEN
          CALL QUIT('Insufficient work space in CC_CMAT. (2b)')
        END IF
       
* read single excitation part of T^A vector:
        IOPT = 1
        CALL CC_RDRSP(LISTXA,IDLSTA,ISYMTA,IOPT,MODEL,
     &                WORK(KT1AMPA),WORK(KDUM))

* read single excitation part of T^B vector:
        IOPT = 1
        CALL CC_RDRSP(LISTXB,IDLSTB,ISYMTB,IOPT,MODEL,
     &                WORK(KT1AMPB),WORK(KDUM))


*=====================================================================*
* A term:
*=====================================================================*

* calculate Gamma^AB: intermediate:
        ISYX4O = MULD2H(ISYOVOV,ISYMAB)

        KX4O    = KEND2
        KXIAJB  = KX4O    + NGAMMA(ISYX4O)
        KENDA3  = KXIAJB  + NT2AM(ISYOVOV)
        LENDA3  = LWORK   - KENDA3

        IF (LENDA3 .LE. 0) THEN
          CALL QUIT('Insufficient work space in CC_CMAT. (A)')
        END IF

* read (ia|jb) integrals from file:
        Call CCG_RDIAJB (WORK(KXIAJB),NT2AM(ISYOVOV))

* calculate double one-index transformed (ik|jl) integrals:
        IOPT = 1
        CALL CCG_4O(WORK(KX4O),ISYX4O,WORK(KXIAJB),ISYOVOV,
     &              WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB,
     &              WORK(KENDA3),LENDA3,IOPT)

* read double excitation part of T^C vector and unpack to full matrix:
        KT2AMPC = KX4O    + NGAMMA(ISYX4O)
        KENDA3  = KT2AMPC + NT2SQ(ISYMTC)
        LENDA3  = LWORK   - KENDA3

        IF (LENDA3 .LE. NT2AM(ISYMTC)) THEN
          CALL QUIT('Insufficient work space in CC_CMAT. (A2)')
        END IF

        IOPT = 2
        CALL CC_RDRSP(LISTXC,IDLSTC,ISYMTC,IOPT,MODEL,
     &                WORK(KDUM),WORK(KENDA3))

        CALL CCLR_DIASCL(WORK(KENDA3),TWO,ISYMTC)

        CALL CC_T2SQ(WORK(KENDA3),WORK(KT2AMPC),ISYMTC)


* contract T^C with Gamma^AB to A term contribution:
        IOPT = 1
        CALL CCRHS_A(WORK(KTHETA2),WORK(KT2AMPC),WORK(KX4O),
     &               WORK(KENDA3),LENDA3,ISYX4O,ISYMTC,IOPT)

        IF (LOCDBG) THEN
          XNORM=DDOT(NGAMMA(ISYX4O),WORK(KX4O),1,WORK(KX4O),1)
          WRITE (LUPRI,*) 'Norm of X4O intermediate:',XNORM
          WRITE (LUPRI,*) MSGDBG,'THETA after A term:'
          WRITE (LUPRI,*) MSGDBG, 'ITRAN/ICYCLE :',ITRAN,ICYCLE
          WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),
     &         ICTRAN(1,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),
     &         ICTRAN(2,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),
     &         ICTRAN(3,ITRAN)
          XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1)
          WRITE (LUPRI,*) MSGDBG, 'Norm of THETA2:',XNORM
          CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1)
          CALL FLSHFO(LUPRI)
        END IF


*=====================================================================*
* B term:
*=====================================================================*
        ISYMKC  = MULD2H(ISYOVOV,ISYMTC)

        KKINTC  = KEND2
        KXIAJB  = KKINTC  + N3ORHF(ISYMKC)
        KT2AMPC = KXIAJB  + NT2SQ(ISYOVOV)
        KENDB3  = KT2AMPC + MAX(NT2AM(ISYMTC),NT2AM(ISYOVOV))
        LENDB3  = LWORK - KENDB3

        IF (LENDB3 .LT. 0) THEN
          CALL QUIT('Insufficient memory in CC_CMAT. (B)')
        END IF

* read (ia|jb) integrals from file and unpack them to a full matrix:
        Call CCG_RDIAJB (WORK(KT2AMPC),NT2AM(ISYOVOV))

        CALL CC_T2SQ(WORK(KT2AMPC),WORK(KXIAJB),ISYOVOV)


* read (ia|jb) integrals from file:
        IOPT = 2
        CALL CC_RDRSP(LISTXC,IDLSTC,ISYMTC,IOPT,MODEL,
     &                WORK(KDUM),WORK(KT2AMPC))
        CALL CCLR_DIASCL(WORK(KT2AMPC),TWO,ISYMTC)

* construct PHI^C intermediate:
        CALL CC_MI(WORK(KKINTC),WORK(KXIAJB),ISYOVOV,
     &             WORK(KT2AMPC),ISYMTC,WORK(KENDB3),LENDB3)

* for CCSD(R12) add correction term:
        IF (CCR12.AND..NOT.(CCS.OR.CC2)) THEN
          KTR12AM = KXIAJB
          KVINT   = KTR12AM + NTR12AM(ISYMTC)
          KENDB3  = KVINT   + NTR12AM(1)
          LENDB3  = LWORK - KENDB3
          IF (LENDB3.LT.0) THEN
            CALL QUIT('Insufficient work space in CC_CMAT (B-R12)')
          END IF

          IOPT = 32
          CALL CC_RDRSP(LISTC,IDLSTC,ISYMTC,IOPT,MODEL,WORK(KDUM),
     &                  WORK(KTR12AM))

          LUNIT = -1
          CALL GPOPEN(LUNIT,FCCR12V,'OLD',' ','UNFORMATTED',
     &                KDUM,.FALSE.)
 9999     READ(LUNIT) IAN
          READ(LUNIT) (WORK(KVINT-1+I), I=1, NTR12AM(1))
          IF (IAN.NE.IANR12) GOTO 9999
          CALL GPCLOSE(LUNIT,'KEEP')

          IOPT = 2
          CALL CC_R12CV(WORK(KKINTC),WORK(KTR12AM),ISYMTC,WORK(KVINT),
     &                  1,IOPT,WORK(KENDB3),LENDB3)

        END IF

* double oneindex-transform PHI^C intermediate to result vector:
        CALL CCC_OVOV(WORK(KTHETA2),ISYRES,WORK(KKINTC),ISYMKC,
     &                WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB,
     &                WORK(KENDB3),LENDB3)



        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG,'THETA after B term:'
          WRITE (LUPRI,*) MSGDBG, 'ITRAN/ICYCLE :',ITRAN,ICYCLE
          WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),
     &         ICTRAN(1,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),
     &         ICTRAN(2,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),
     &         ICTRAN(3,ITRAN)
          XNORM=DDOT(N3ORHF(ISYMKC),WORK(KKINTC),1,WORK(KKINTC),1)
          WRITE (LUPRI,*) MSGDBG, 'Norm of KINTC:',XNORM
          XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1)
          WRITE (LUPRI,*) MSGDBG, 'Norm of THETA2:',XNORM
          CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1)
          CALL FLSHFO(LUPRI)
        END IF

*=====================================================================*
* C term: requires 5 x 1/2 V^2 O^2  !!!
* D term: requires 5 x 1/2 V^2 O^2  !!!
*=====================================================================*
        ISYCBAR = MULD2H(ISYMTC,ISYOVOV)
        ISYDBAR = MULD2H(ISYMTC,ISYOVOV)

        KXIAJB  = KEND2
        KT2AMPC = KXIAJB  + NT2SQ(ISYOVOV)
        KCBAR   = KT2AMPC + MAX(NT2AM(ISYMTC),NT2AM(ISYOVOV))
        KEND3   = KCBAR   + NT2SQ(ISYCBAR)
        LEND3   = LWORK - KEND3

        KDBAR   = KCBAR

        IF (LEND3 .LT. 0) THEN
          CALL QUIT('Insufficient work space in CC_BMAT. (C)')
        END IF

* read (ia|jb) and square them:
        Call CCG_RDIAJB (WORK(KT2AMPC),NT2AM(ISYOVOV))

        Call CC_T2SQ (WORK(KT2AMPC), WORK(KXIAJB), ISYOVOV)

* read double excitation part of T^C vector:
        IOPT = 2
        CALL CC_RDRSP(LISTXC,IDLSTC,ISYMTC,IOPT,MODEL,
     &                WORK(KDUM),WORK(KT2AMPC)          )

        CALL CCLR_DIASCL(WORK(KT2AMPC),TWO,ISYMTC)

* calculate CBAR^C intermediate:
        IOPT   = 1
        CBAFIL = '--------'
        LUCBAR = -1
        IOFFCD = -1
        CALL CCB_CDBAR('C',WORK(KXIAJB),ISYOVOV, WORK(KT2AMPC),ISYMTC,
     &                     WORK(KCBAR), ISYCBAR, WORK(KEND3),LEND3,
     &                     CBAFIL, LUCBAR, IOFFCD, IOPT)

* double transform CBAR^C intermediate to C term of the result vector:
        CALL CCB_22CD(WORK(KTHETA2),ISYRES,WORK(KCBAR),ISYCBAR, 
     &                WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB,
     &                'C', WORK(KEND3),LEND3)


* calculate DBAR^D intermediate:
        IOPT   = 1
        DBAFIL = '--------'
        LUDBAR = -1
        IOFFCD = -1
        CALL CCB_CDBAR('D',WORK(KXIAJB),ISYOVOV, WORK(KT2AMPC),ISYMTC,
     &                     WORK(KDBAR), ISYDBAR, WORK(KEND3),LEND3,
     &                     DBAFIL, LUDBAR, IOFFCD, IOPT)

* double transform DBAR^C intermediate to D term of the result vector:
        CALL CCB_22CD(WORK(KTHETA2),ISYRES,WORK(KDBAR),ISYDBAR, 
     &                WORK(KT1AMPA),ISYMTA,WORK(KT1AMPB),ISYMTB,
     &                'D', WORK(KEND3),LEND3)


        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG,'THETA after C and D terms:'
          WRITE (LUPRI,*) MSGDBG, 'ITRAN/ICYCLE :',ITRAN,ICYCLE
          WRITE (LUPRI,*) MSGDBG, 'LISTA/IDLSTA:',LISTA(1:2),
     &         ICTRAN(1,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTB/IDLSTB:',LISTB(1:2),
     &         ICTRAN(2,ITRAN)
          WRITE (LUPRI,*) MSGDBG, 'LISTC/IDLSTC:',LISTC(1:2),
     &         ICTRAN(3,ITRAN)
          XNORM=DDOT(NT2AM(ISYRES),WORK(KTHETA2),1,WORK(KTHETA2),1)
          WRITE (LUPRI,*) MSGDBG, 'Norm of THETA2:',XNORM
          CALL CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYRES,1,1)
          CALL FLSHFO(LUPRI)
        END IF

*---------------------------------------------------------------------*
* End loop over cyclic permutations, save result vector
*---------------------------------------------------------------------*
        END DO ! ICYCLE

        KTHETA0 = -999999

        IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN
          CALL PUTWA2(LUCMAT,FILCMA,WORK(KTHETA1),
     &                ICTRAN(4,ITRAN),NT1AM(ISYRES)+NT2AM(ISYRES))
c       ELSE IF (IOPTRES.EQ.3) THEN
c         CALL CC_WRRSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW,
c    &                  WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
c    &                  WORK(KEND1),LEND1)
        ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN
          CALL CC_WARSP(FILCMA,ICTRAN(4,ITRAN),ISYRES,IOPTW,MODELW,
     &                  WORK(KTHETA0),WORK(KTHETA1),WORK(KTHETA2),
     &                  WORK(KEND1),LEND1)
        ELSE IF (IOPTRES.EQ.5) THEN
          CALL CCDOTRSP(ICDOTS,CCONS,IOPTW,FILCMA,ITRAN,NCTRAN,MXVEC,
     &                  WORK(KTHETA1),WORK(KTHETA2),ISYRES,
     &                  WORK(KEND1),LEND1)
        ELSE
          CALL QUIT('Error in CC_CMAT.')
        END IF
*---------------------------------------------------------------------*
* End loop over C matrix transformations
*---------------------------------------------------------------------*
      END DO

        IF (LOCDBG) THEN
          WRITE (LUPRI,*) MSGDBG,'N^6 term section finished...'
          CALL FLSHFO(LUPRI)
        END IF

      END IF ! (.NOT. (CCS .OR. CC2))
*=====================================================================*
* End of the N^6 term section
*=====================================================================*

*=====================================================================*
* restore result vectors and clean up and return:
*=====================================================================*
8888   CONTINUE

*---------------------------------------------------------------------*
* if IOPTRES=1 and enough work space available, read result
* vectors back into memory:
*---------------------------------------------------------------------*

* check size of work space:
      IF (IOPTRES .EQ. 1) THEN
        LENALL = IADRTH-1
        IF (LENALL .GT. LWORK) IOPTRES = 0
      END IF

* read the result vectors back into memory:
      IF (IOPTRES .EQ. 1) THEN

        CALL GETWA2(LUCMAT,FILCMA,WORK(1),1,LENALL)

        IF (LOCDBG) THEN
          DO ITRAN = 1, NCTRAN
            IF (ITRAN.LT.NCTRAN) THEN
              LEN     = ICTRAN(4,ITRAN+1)-ICTRAN(4,ITRAN)
            ELSE
              LEN     = IADRTH-ICTRAN(4,NCTRAN)
            END IF
            KTHETA1 = ICTRAN(4,ITRAN)
            XNORM   = DDOT(LEN, WORK(KTHETA1),1, WORK(KTHETA1),1)
            WRITE (LUPRI,*) 'Read C matrix transformation nb. ',ITRAN
            WRITE (LUPRI,*) 'Adress, length, NORM:',ICTRAN(4,NCTRAN),
     &           LEN,XNORM
          END DO
          CALL FLSHFO(LUPRI)
        END IF
      END IF 

*---------------------------------------------------------------------*
* close C matrix file & return
*---------------------------------------------------------------------*
* check return option for the result vectors:
      IF (IOPTRES .EQ. 0 .OR. IOPTRES .EQ. 1) THEN
        CALL WCLOSE2(LUCMAT,FILCMA,'KEEP')

      ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4 .OR. IOPTRES.EQ.5) THEN
        CONTINUE
      ELSE
        CALL QUIT('Illegal value of IOPTRES in CC_CMAT.')
      END IF


*=====================================================================*

      CALL QEXIT('CC_CMAT')

      RETURN
      END
*=====================================================================*
*            END OF SUBROUTINE CC_CMAT
*=====================================================================*

*---------------------------------------------------------------------*
c/* Deck CC_CTST */
*=====================================================================*
       SUBROUTINE CC_CTST(WORK,LWORK)
#if defined (IMPLICIT_NONE)
      IMPLICIT NONE  
#else
#  include "implicit.h"
#endif
#include "priunit.h"
#include "ccsdinp.h"
#include "ccsdsym.h"
#include "ccorb.h"

* local parameters:
      CHARACTER MSGDBG*(18)
      PARAMETER (MSGDBG='[debug] CC_CTST> ')

      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.)
      INTEGER MXCTRAN
      PARAMETER (MXCTRAN = 2)

      INTEGER LWORK
      DOUBLE PRECISION WORK(LWORK) 
      DOUBLE PRECISION DDOT, RDUM 

      CHARACTER*(3) LISTA, LISTB, LISTC, LISTD
      CHARACTER*(8) FILCMA, LABELA
      CHARACTER*(10) MODEL
      INTEGER IOPTRES, IDUM
      INTEGER ICTRAN(4,MXCTRAN), NCTRAN
      INTEGER IDLSTA, IDLSTB, IDLSTC, ISYMA, ISYMB, ISYMC, ISYMABC
      INTEGER KTHETA1, KTHETA2, KT1AMPC, KT2AMPC, KRESLT1, KRESLT2
      INTEGER KEND1, LEND1, IOPT, ISYMD, KT1AMPD, KT2AMPD, IDLSTD

* external function:
      INTEGER IR1TAMP
      INTEGER IL1ZETA
      INTEGER ILSTSYM


      CALL QENTER('CC_CTST')


*---------------------------------------------------------------------*
* call C matrix transformation:
*---------------------------------------------------------------------*
      LISTA = 'R1'
      LISTB = 'R1'
      LISTC = 'R1'
      LISTD = 'L1'
      IDLSTA = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1)
      IDLSTB = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1)
      IDLSTC = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,1)
      IDLSTD = IL1ZETA('ZDIPLEN ',.FALSE.,0.0D0,1)
      ICTRAN(1,1) = IDLSTA
      ICTRAN(2,1) = IDLSTB
      ICTRAN(3,1) = IDLSTC
      NCTRAN = 1

      IOPTRES = 1
      FILCMA  = 'CCCMAT'

      WRITE (LUPRI,*) 'CCTST: call CC_CMAT:'

      CALL CC_CMAT(ICTRAN,  NCTRAN,  LISTA,  LISTB, LISTC,
     &             IOPTRES, FILCMA, IDUM, RDUM, 0, WORK, LWORK )


      ISYMA  = ILSTSYM(LISTA,IDLSTA)
      ISYMB  = ILSTSYM(LISTB,IDLSTB)
      ISYMC  = ILSTSYM(LISTC,IDLSTC)
      ISYMD  = ILSTSYM(LISTD,IDLSTD)
      ISYMABC = MULD2H(MULD2H(ISYMA,ISYMB),ISYMC)

      KTHETA1 = ICTRAN(4,1)
      KTHETA2 = KTHETA1 + NT1AM(ISYMABC)
      KEND1   = KTHETA2 + NT2AM(ISYMABC)
      LEND1   = LWORK - KEND1

      IF (NSYM.EQ.1 .AND. LOCDBG) THEN
        KT1AMPC = KTHETA2 + NT2AM(ISYMABC)
        KT2AMPC = KT1AMPC + NT1AM(ISYMC)
        KRESLT1 = KT2AMPC + NT2AM(ISYMC)
        KRESLT2 = KRESLT1 + NT1AM(ISYMABC)
        KEND1   = KRESLT2 + NT2AM(ISYMABC)
        LEND1   = LWORK - KEND1

        IF (LEND1 .LT. 0) THEN
          CALL QUIT('Insufficient work space in CC_CTST.')
        END IF

        IOPT = 3
        Call CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL,
     &                WORK(KT1AMPC),WORK(KT2AMPC))

        ! zero singles or doubles C vector:
C       CALL DZERO(WORK(KT1AMPC),NT1AM(ISYMC))
C       CALL DZERO(WORK(KT2AMPC),NT2AM(ISYMC))
        CALL DZERO(WORK(KRESLT1),NT1AM(ISYMABC)+NT2AM(ISYMABC))
        IPRINT  = 5

        CALL CC_FDC(NT1AM(ISYMABC),NT2AM(ISYMABC),
     >              LISTA,IDLSTA,LISTB,IDLSTB,  
     >              WORK(KT1AMPC), WORK(KRESLT1),
     >              WORK(KEND1), LEND1)

        IPRINT  = 0

        IF (.TRUE.) THEN
          WRITE (LUPRI,*) 
     *          'LISTA, IDLSTA, ISYMA:',LISTA(1:2),IDLSTA,ISYMA
          WRITE (LUPRI,*) 
     *          'LISTB, IDLSTB, ISYMB:',LISTB(1:2),IDLSTB,ISYMB
          WRITE (LUPRI,*) 
     *          'LISTC, IDLSTC, ISYMC:',LISTC(1:2),IDLSTC,ISYMC
          WRITE (LUPRI,*) 'ISYMABC:',ISYMABC
          WRITE (LUPRI,*)
          WRITE (LUPRI,*) 'finite difference Theta vector:'
          Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYMABC,1,1)
          WRITE (LUPRI,*) 'analytical Theta vector:'
          Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYMABC,1,1)
        END IF

        KT1AMPD = KEND1
        KT2AMPD = KT1AMPD + NT1AM(ISYMABC)
        KEND1   = KT2AMPD + NT2AM(ISYMABC)
        LEND1   = LWORK - KEND1
        IF (LEND1 .LT. 0) THEN
          CALL QUIT('Insufficient work space in CC_CTST.')
        END IF

        IOPT = 3
        Call CC_RDRSP(LISTD,IDLSTD,ISYMD,IOPT,MODEL,
     &                WORK(KT1AMPD),WORK(KT2AMPD))
        CALL CCLR_DIASCL(WORK(KT2AMPD),0.5d0,ISYMD)

        WRITE (LUPRI,*) 'Dot product with T^D for analytical C matrix:',
     >   DDOT(NT1AM(ISYMD),WORK(KT1AMPD),1,WORK(KTHETA1),1)+
     >   DDOT(NT2AM(ISYMD),WORK(KT2AMPD),1,WORK(KTHETA2),1) 

        WRITE (LUPRI,*) 'Dot product with T^D for numerical C matrix:',
     >   DDOT(NT1AM(ISYMD),WORK(KT1AMPD),1,WORK(KRESLT1),1)+
     >   DDOT(NT2AM(ISYMD),WORK(KT2AMPD),1,WORK(KRESLT2),1) 



        Call DAXPY(NT1AM(ISYMABC),-1.0d0,WORK(KTHETA1),1,
     &                                  WORK(KRESLT1),1)
        IF (.NOT.CCS) THEN
          Call DAXPY(NT2AM(ISYMABC),-1.0d0,WORK(KTHETA2),1,
     &                                    WORK(KRESLT2),1)
        ELSE
          Call DZERO(WORK(KRESLT2),NT2AM(ISYMABC))
        END IF

        WRITE (LUPRI,*) 'Norm of difference between analytical THETA '
     >           // 'vector and the numerical result:'
        WRITE (LUPRI,*) 'singles excitation part:',
     >   DSQRT(DDOT(NT1AM(ISYMA),WORK(KRESLT1),1,WORK(KRESLT1),1))
        WRITE (LUPRI,*) 'double excitation part: ',
     >   DSQRT(DDOT(NT2AM(ISYMA),WORK(KRESLT2),1,WORK(KRESLT2),1))

        WRITE (LUPRI,*) 'difference vector:'
        Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYMABC,1,1)

        CALL FLSHFO(LUPRI)


      ELSE IF (NSYM.NE.1 .AND. LOCDBG) THEN
        WRITE (LUPRI,*) 'analytical Theta vector:'
        Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYMABC,1,1)

        WRITE (LUPRI,*) 'CC_CTST> can not calculate finite '//
     &       'difference C matrix'
        WRITE (LUPRI,*) 'CC_CTST> with symmetry.'

        KT1AMPD = KEND1
        KT2AMPD = KT1AMPD + NT1AM(ISYMABC)
        KEND1   = KT2AMPD + NT2AM(ISYMABC)
        LEND1   = LWORK - KEND1
        IF (LEND1 .LT. 0) THEN
          CALL QUIT('Insufficient work space in CC_CTST.')
        END IF

        IOPT = 3
        Call CC_RDRSP(LISTD,IDLSTD,ISYMD,IOPT,MODEL,
     &                WORK(KT1AMPD),WORK(KT2AMPD))
        CALL CCLR_DIASCL(WORK(KT2AMPD),0.5d0,ISYMD)

        WRITE (LUPRI,*) 'Dot product with T^D for analytical C matrix:',
     >   DDOT(NT1AM(ISYMD),WORK(KT1AMPD),1,WORK(KTHETA1),1)+
     >   DDOT(NT2AM(ISYMD),WORK(KT2AMPD),1,WORK(KTHETA2),1) 

      END IF

      CALL QEXIT('CC_CTST')

      RETURN
      END 
*=====================================================================*
*---------------------------------------------------------------------*
c/*  Deck CCC_OVOV */
*=====================================================================*
      SUBROUTINE CCC_OVOV(THETA2, ISYRES, XKINT,  ISYXKI, 
     &                    T1AMPA, ISYMTA, T1AMPB, ISYMTB,
     &                    WORK, LWORK )
*---------------------------------------------------------------------*
*
*     Purpose: double one-index transform XKINT_{iljk} intermediate 
*              with T^A and T^B single-excitation amplitudes 
*
*              THETA2(ai,bj) += P^{AB}  t^A_{ak} t^B_{bl} XKINT_{iljk}
*
*              needed for CCSD part of C matrix transformation
*
*     Christof Haettig, Maj 1997
*=====================================================================*
#if defined (IMPLICIT_NONE)
      IMPLICIT NONE
#else
# include "implicit.h"
#endif
 
#include "priunit.h"
#include "ccsdsym.h"
#include "ccorb.h"
 
      DOUBLE PRECISION ZERO, ONE, TWO
      PARAMETER (ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0)

      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.)

      INTEGER ISYMTA, ISYMTB, ISYXKI, ISYRES
      INTEGER LWORK

      DOUBLE PRECISION THETA2(*)     ! dimension (NT2AM(ISYRES))
      DOUBLE PRECISION XKINT(*)      ! dimension (N3ORHF(ISYXKI))
      DOUBLE PRECISION T1AMPA(*)     ! dimension (NT1AM(ISYMTA))
      DOUBLE PRECISION T1AMPB(*)     ! dimension (NT1AM(ISYMTB))
      DOUBLE PRECISION WORK(LWORK)
      DOUBLE PRECISION SUM1, DDOT, SUM2

      INTEGER ISYMA, ISYMI, ISYMB, ISYMJ, ISYMK, ISYML
      INTEGER ISYMAI, ISYMBJ, ISYMJL, ISYJLI, NAIBJ
      INTEGER KSCR1, KSCR2, KEND1, LEND1, KEND2, LEND2, KOFF1, KOFF2
      INTEGER KOFFX, KOFFT, NTOTJLI, NTOTA, NTOTB, NTOTJ, NAI, NBJ


* statement function:
      INTEGER INDEX
      INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J

      CALL QENTER('CCC_OVOV')

*---------------------------------------------------------------------*
* begin:
*---------------------------------------------------------------------*
* check symmetries:
      IF (MULD2H(ISYMTA,ISYMTB) .NE. MULD2H(ISYRES,ISYXKI)) THEN
        CALL QUIT('Symmetry mismatch in CCC_OVOV.')
      END IF

      IF (LOCDBG) THEN
        WRITE (LUPRI,*) 'CCC_OVOV> ISYRES:', ISYRES
        WRITE (LUPRI,*) 'CCC_OVOV> ISYXKI:', ISYXKI
        WRITE (LUPRI,*) 'CCC_OVOV> ISYMTA:', ISYMTA
        WRITE (LUPRI,*) 'CCC_OVOV> ISYMTB:', ISYMTB
        WRITE (LUPRI,*) 'CCC_OVOV> XKINT:'
C       WRITE (LUPRI,'(5F14.6)') XKINT
      END IF

      SUM1 = 0.0d0

*---------------------------------------------------------------------*
* loop over A indeces:
*---------------------------------------------------------------------*
      DO ISYMA = 1, NSYM
        ISYMK  = MULD2H(ISYMA,ISYMTA)
        ISYJLI = MULD2H(ISYXKI,ISYMK)

      IF (NRHF(ISYMK).NE.0) THEN


        KSCR1 = 1
        KEND1 = KSCR1 + NMAIJK(ISYJLI)
        LEND1 = LWORK - KEND1
    
        IF (LEND1.LT.0) THEN
          CALL QUIT('Insufficient memory in CCC_OVOV.')
        END IF

      DO A = 1, NVIR(ISYMA)

*---------------------------------------------------------------------*
* transform K index:  sum_k XKINT_{jli;k} * t^A(ka)
*---------------------------------------------------------------------*
          KOFFX = I3ORHF(ISYJLI,ISYMK) + 1
          KOFFT = IT1AM(ISYMA,ISYMK) + A

          NTOTJLI = MAX(1,NMAIJK(ISYJLI))
          NTOTA   = MAX(1,NVIR(ISYMA))

          Call DGEMV('N', NMAIJK(ISYJLI), NRHF(ISYMK),
     &                ONE,  XKINT(KOFFX),  NTOTJLI, 
     &                      T1AMPA(KOFFT), NTOTA,
     &                ZERO, WORK(KSCR1),   1          )

*---------------------------------------------------------------------*
* loop over I indeces:
*---------------------------------------------------------------------*
      DO ISYMI = 1, NSYM
        ISYMJL = MULD2H(ISYJLI,ISYMI)
        ISYMAI = MULD2H(ISYMA,ISYMI)
        ISYMBJ = MULD2H(ISYRES,ISYMAI)

        KSCR2 = KEND1
        KEND2 = KSCR2 + NT1AM(ISYMBJ)
        LEND2 = LWORK - KEND2
   
        IF (LEND2.LT.0) THEN
          CALL QUIT('Insufficient memory in CCC_OVOV.')
        END IF


      DO I = 1, NRHF(ISYMI)

*---------------------------------------------------------------------*
* transform L index:  sum_l t^B(bl) * SCR_{jl;i}
*---------------------------------------------------------------------*
      DO ISYMB = 1, NSYM
        ISYMJ = MULD2H(ISYMBJ,ISYMB)
        ISYML = MULD2H(ISYMTB,ISYMB)
        IF (MULD2H(ISYMJ,ISYML).NE.ISYMJL)
     &       CALL QUIT('Error in CCC_OVOV (1).')
        
        KOFF1 = KSCR1 + IMAIJK(ISYMJL,ISYMI) +
     &            NMATIJ(ISYMJL)*(I-1) + IMATIJ(ISYMJ,ISYML)
        KOFF2 = KSCR2 + IT1AM(ISYMB,ISYMJ)
        KOFFT = IT1AM(ISYMB,ISYML) + 1

        NTOTJ = MAX(1,NRHF(ISYMJ))
        NTOTB = MAX(1,NVIR(ISYMB))

        Call DGEMM('N','T',NVIR(ISYMB),NRHF(ISYMJ),NRHF(ISYML),
     &              ONE, T1AMPB(KOFFT), NTOTB, WORK(KOFF1), NTOTJ,
     &             ZERO, WORK(KOFF2),   NTOTB )
      END DO

*---------------------------------------------------------------------*
* store result in THETA2 vector:
*---------------------------------------------------------------------*
      NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A

      IF      (ISYMAI .EQ. ISYMBJ) THEN

        WORK(KSCR2-1+NAI) = TWO * WORK(KSCR2-1+NAI)

        DO NBJ = 1, NT1AM(ISYMBJ)
          NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NAI,NBJ)
          THETA2(NAIBJ) = THETA2(NAIBJ) + WORK(KSCR2-1+NBJ)
        END DO

      ELSE IF (ISYMAI .LT. ISYMBJ) THEN

        DO NBJ = 1, NT1AM(ISYMBJ)
          NAIBJ = IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1)+NAI
          THETA2(NAIBJ) = THETA2(NAIBJ) + WORK(KSCR2-1+NBJ)
        END DO

      ELSE IF (ISYMBJ .LT. ISYMAI) THEN

        DO NBJ = 1, NT1AM(ISYMBJ)
          NAIBJ = IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMBJ)*(NAI-1)+NBJ
          THETA2(NAIBJ) = THETA2(NAIBJ) + WORK(KSCR2-1+NBJ)
        END DO

      ELSE
        CALL QUIT('Error in CCC_OVOV (2).')
      END IF

      

*---------------------------------------------------------------------*
* end loops over I and A:
*---------------------------------------------------------------------*
      END DO ! I
      END DO ! ISYMI

      END DO ! A
      END IF ! NRHF(ISYMK).NE.0
      END DO ! ISYMA

      CALL QEXIT('CCC_OVOV')

      RETURN
      END
*---------------------------------------------------------------------*
*                   END OF SUBROUTINE CCC_OVOV                        *
*---------------------------------------------------------------------*
